Signal Decay¶
import os
import imageio
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from IPython.display import Image
from nilearn.glm import first_level
from scipy import signal
from book_utils import predict_bold_signal
sns.set_style("whitegrid")
plt.rcParams.update({
"text.usetex": True,
"font.family": "sans-serif",
"font.sans-serif": ["Helvetica"]
})
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/nilearn/datasets/__init__.py:96: FutureWarning: Fetchers from the nilearn.datasets module will be updated in version 0.9 to return python strings instead of bytes and Pandas dataframes instead of Numpy arrays.
"Numpy arrays.", FutureWarning)
/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/site-packages/nilearn/glm/__init__.py:56: FutureWarning: The nilearn.glm module is experimental. It may change in any future release of Nilearn.
'It may change in any future release of Nilearn.', FutureWarning)
import plotly.graph_objects as go
import numpy as np
# Create figure
fig = go.Figure()
# Add traces, one for each slider step
for step in np.arange(0, 5, 0.1):
fig.add_trace(
go.Scatter(
visible=False,
line=dict(color="#00CED1", width=6),
name="𝜈 = " + str(step),
x=np.arange(0, 10, 0.01),
y=np.sin(step * np.arange(0, 10, 0.01))))
# Make 10th trace visible
fig.data[10].visible = True
# Create and add slider
steps = []
for i in range(len(fig.data)):
step = dict(
method="update",
args=[{"visible": [False] * len(fig.data)},
{"title": "Slider switched to step: " + str(i)}], # layout attribute
)
step["args"][0]["visible"][i] = True # Toggle i'th trace to "visible"
steps.append(step)
sliders = [dict(
active=10,
currentvalue={"prefix": "Frequency: "},
pad={"t": 50},
steps=steps
)]
fig.update_layout(
sliders=sliders
)
fig.show()
Signal decays as echo time increases¶
data_dir = os.path.abspath("../data/")
files = sorted(glob(os.path.join(data_dir, "sub-pilot_task-checkerboard_run-2_*_bold.nii.gz")))
echo_times = np.array([9.58, 21.95, 34.32, 46.69, 59.06, 71.43, 83.8, 96.17])
img = image.index_img(files[0], 0)
data = img.get_fdata()
vmax = np.max(data)
idx = np.vstack(np.where(data > 1000))
min_x, min_y, min_z = np.min(idx, axis=1)
max_x, max_y, max_z = np.max(idx, axis=1)
imgs = []
for f in files:
img = image.index_img(f, 0)
data = img.get_fdata()
data = data[min_x:max_x, min_y:max_y, min_z:max_z]
img = nib.Nifti1Image(data, img.affine, img.header)
imgs.append(img)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/tmp/ipykernel_2056/1903662462.py in <module>
1 data_dir = os.path.abspath("../data/")
2
----> 3 files = sorted(glob(os.path.join(data_dir, "sub-pilot_task-checkerboard_run-2_*_bold.nii.gz")))
4 echo_times = np.array([9.58, 21.95, 34.32, 46.69, 59.06, 71.43, 83.8, 96.17])
5
NameError: name 'glob' is not defined
plt.style.use("dark_background")
fig, axes = plt.subplots(figsize=(26, 4), ncols=len(files))
for i_echo, img in enumerate(imgs):
te = echo_times[i_echo]
if i_echo == 0:
data = img.get_fdata()
vmax = np.max(data)
plotting.plot_epi(
img,
cut_coords=[0],
display_mode="z",
annotate=False,
draw_cross=False,
axes=axes[i_echo],
figure=fig,
vmax=vmax,
cmap="gray",
)
axes[i_echo].set_title(f"TE={te}ms", fontsize=20, pad=0)
fig.show()
The impact of \(T_{2}^{*}\) and \(S_{0}\) fluctuations on BOLD signal¶
# Constants
OUT_DIR = "figures/"
os.makedirs(OUT_DIR, exist_ok=True)
# Simulate data
MULTIECHO_TES = np.array([15, 30, 45, 60, 75, 90])
SINGLEECHO_TE = np.array([30])
# For a nice, smooth curve
FULLCURVE_TES = np.arange(0, 101, 1)
n_echoes = len(FULLCURVE_TES)
pal = sns.color_palette('cubehelix', 8)
Simulate \(T_{2}^{*}\) and \(S_{0}\) fluctuations¶
# Simulate data
# We'll convolve with HRF just for smoothness
hrf = first_level.spm_hrf(1, oversampling=1)
N_VOLS = 21
SCALING_FRACTION = 0.1 # used to scale standard deviation
MEAN_T2S = 35
MEAN_S0 = 16000
# simulate the T2*/S0 time series
# The original time series will be a random time series from a normal distribution,
# convolved with the HRF
ts = np.random.normal(loc=0, scale=1, size=(N_VOLS+20,))
ts = signal.convolve(ts, hrf)[20:N_VOLS+20]
ts *= SCALING_FRACTION / np.std(ts)
ts -= np.mean(ts)
t2s_ts = (ts * MEAN_T2S) + MEAN_T2S
s0_ts = (ts * MEAN_S0) + MEAN_S0
Plot BOLD signal decay for a standard single-echo scan¶
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [30])
fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
FULLCURVE_TES,
fullcurve_signal[:, 0],
alpha=0.5,
color="black",
)
ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.show()
Plot BOLD signal decay and BOLD contrast¶
fullcurve_signal_active = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [40])
fullcurve_signal_inactive = predict_bold_signal(FULLCURVE_TES, [MEAN_S0], [20])
fig, ax = plt.subplots(figsize=(14, 7.5))
ax.plot(
FULLCURVE_TES,
fullcurve_signal_active[:, 0],
alpha=0.5,
color="red",
label="High Activity",
)
ax.plot(
FULLCURVE_TES,
fullcurve_signal_inactive[:, 0],
alpha=0.5,
color="blue",
label="Low Activity",
)
ax.set_ylabel("Recorded BOLD signal", fontsize=24)
ax.set_xlabel("Echo Time (ms)", fontsize=24)
ax.set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
ax.set_xlim(0, np.max(FULLCURVE_TES))
ax.legend(fontsize=20)
ax.tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.show()
Plot single-echo data resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how fMRI data fluctuates over time.
Plot single-echo data and the curve resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how single-echo data is a sample from a signal decay curve.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
out_file = os.path.join(OUT_DIR, "fluctuations_single-echo_with_curve.gif")
if os.path.isfile(out_file):
os.remove(out_file)
filenames = []
for i_vol in range(N_VOLS):
filename = f"se_{i_vol}.png"
fig, axes = plt.subplots(
nrows=2,
figsize=(14, 10),
gridspec_kw={"height_ratios": [1, 3]}
)
axes[0].plot(
singleecho_signal[0, :],
color="black",
zorder=0
)
axes[0].scatter(
i_vol,
singleecho_signal[:, i_vol],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].plot(
FULLCURVE_TES,
fullcurve_signal[:, i_vol],
alpha=0.5,
color="black",
zorder=0,
)
axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, i_vol],
color="orange",
s=150,
alpha=1.,
label="Single-Echo Signal",
zorder=1,
)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.savefig(filename)
plt.close(fig)
filenames.append(filename)
# build gif
with imageio.get_writer(out_file, mode="I") as writer:
for filename in filenames:
image = imageio.imread(filename)
writer.append_data(image)
# Remove files
for filename in filenames:
os.remove(filename)
with open(out_file, "rb") as file:
display(Image(file.read(), width=600))
Plot single-echo data, the curve, and the \(S_{0}\) and \(T_{2}^{*}\) values resulting from \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how changes in fMRI data can be driven by both S0 and T2* fluctuations.
fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, t2s_ts)
singleecho_signal = fullcurve_signal[SINGLEECHO_TE, :]
out_file = os.path.join(OUT_DIR, "fluctuations_single-echo_with_curve_and_t2s_s0.gif")
if os.path.isfile(out_file):
os.remove(out_file)
filenames = []
for i_vol in range(N_VOLS):
filename = f"se_{i_vol}.png"
fig, axes = plt.subplots(
nrows=2,
figsize=(14, 10),
gridspec_kw={"height_ratios": [1, 3]}
)
axes[0].plot(
singleecho_signal[0, :],
color="black",
zorder=0,
)
axes[0].scatter(
i_vol,
singleecho_signal[:, i_vol],
color="orange",
s=150,
label="Single-Echo Signal",
zorder=1,
)
axes[0].set_ylabel("Signal", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].scatter(
SINGLEECHO_TE,
singleecho_signal[:, i_vol],
color="orange",
s=150,
alpha=1.,
label="Single-Echo Signal",
zorder=3,
)
axes[1].plot(
FULLCURVE_TES,
fullcurve_signal[:, i_vol],
alpha=0.5,
color="black",
zorder=0,
)
t2s_value = predict_bold_signal(
np.array([t2s_ts[i_vol]]),
np.array([s0_ts[i_vol]]),
np.array([t2s_ts[i_vol]])
)[0]
axes[1].scatter(
t2s_ts[i_vol],
t2s_value,
marker="*",
color="blue",
s=500,
alpha=0.5,
label="$T_{2}^{*}$",
zorder=1,
)
axes[1].scatter(
0,
s0_ts[i_vol],
marker="*",
color="red",
s=500,
alpha=0.5,
label="$S_{0}$",
clip_on=False,
zorder=2,
)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.savefig(filename)
plt.close(fig)
filenames.append(filename)
# build gif
with imageio.get_writer(out_file, mode="I") as writer:
for filename in filenames:
image = imageio.imread(filename)
writer.append_data(image)
# Remove files
for filename in filenames:
os.remove(filename)
with open(out_file, "rb") as file:
display(Image(file.read(), width=600))
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations¶
This shows how fluctuations in S0 and T2* produce different patterns in the full signal decay curves.
s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
out_file = os.path.join(OUT_DIR, "fluctuations_t2s_s0.gif")
if os.path.isfile(out_file):
os.remove(out_file)
filenames = []
for i_vol in range(N_VOLS):
filename = f"total_{i_vol}.png"
fig, axes = plt.subplots(
nrows=2,
figsize=(14, 10),
gridspec_kw={"height_ratios": [1, 3]}
)
axes[0].plot(ts, color="black")
axes[0].scatter(
i_vol,
ts[i_vol],
color="purple",
s=150,
)
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, i_vol],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
)
axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, i_vol],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.savefig(filename)
plt.close(fig)
filenames.append(filename)
# build gif
with imageio.get_writer(out_file, mode="I") as writer:
for filename in filenames:
image = imageio.imread(filename)
writer.append_data(image)
# Remove files
for filename in filenames:
os.remove(filename)
with open(out_file, "rb") as file:
display(Image(file.read(), width=600))
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting single-echo data¶
This shows how single-echo data, on its own, cannot distinguish between S0 and T2* fluctuations.
s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
s0based_singleecho_signal = s0based_fullcurve_signal[SINGLEECHO_TE, :]
t2sbased_singleecho_signal = t2sbased_fullcurve_signal[SINGLEECHO_TE, :]
out_file = os.path.join(OUT_DIR, "fluctuations_t2s_s0_single-echo.gif")
if os.path.isfile(out_file):
os.remove(out_file)
filenames = []
for i_vol in range(N_VOLS):
filename = f"total_{i_vol}.png"
fig, axes = plt.subplots(
nrows=2,
figsize=(14, 10),
gridspec_kw={"height_ratios": [1, 3]}
)
axes[0].plot(
ts,
color="black",
zorder=0,
)
axes[0].scatter(
i_vol,
ts[i_vol],
color="purple",
s=150,
zorder=1,
)
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, i_vol],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
zorder=0,
)
axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, i_vol],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
zorder=1,
)
axes[1].scatter(
SINGLEECHO_TE,
s0based_singleecho_signal[:, i_vol],
color="red",
s=150,
alpha=1.,
label="$S_0$-Driven Single-Echo Signal",
zorder=2,
)
axes[1].scatter(
SINGLEECHO_TE,
t2sbased_singleecho_signal[:, i_vol],
color="blue",
s=150,
alpha=1.,
label="$T_{2}^{*}$-Driven Single-Echo Signal",
zorder=3,
)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.savefig(filename)
plt.close(fig)
filenames.append(filename)
# build gif
with imageio.get_writer(out_file, mode="I") as writer:
for filename in filenames:
image = imageio.imread(filename)
writer.append_data(image)
# Remove files
for filename in filenames:
os.remove(filename)
with open(out_file, "rb") as file:
display(Image(file.read(), width=600))
Plot \(S_{0}\) and \(T_{2}^{*}\) fluctuations and resulting multi-echo data¶
This shows how S0 and T2* fluctuations produce different patterns in multi-echo data.
s0based_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, s0_ts, np.full(N_VOLS, MEAN_T2S))
t2sbased_fullcurve_signal = predict_bold_signal(FULLCURVE_TES, np.full(N_VOLS, MEAN_S0), t2s_ts)
s0based_multiecho_signal = s0based_fullcurve_signal[MULTIECHO_TES, :]
t2sbased_multiecho_signal = t2sbased_fullcurve_signal[MULTIECHO_TES, :]
out_file = os.path.join(OUT_DIR, "fluctuations_t2s_s0_multi-echo.gif")
if os.path.isfile(out_file):
os.remove(out_file)
filenames = []
for i_vol in range(N_VOLS):
filename = f"total_{i_vol}.png"
fig, axes = plt.subplots(
nrows=2,
figsize=(14, 10),
gridspec_kw={"height_ratios": [1, 3]}
)
axes[0].plot(
ts,
color="black",
zorder=0,
)
axes[0].scatter(
i_vol,
ts[i_vol],
color="purple",
s=150,
zorder=1,
)
axes[0].set_ylabel("$S_0$/$T_{2}^{*}$", fontsize=24)
axes[0].set_xlabel("Volume", fontsize=24)
axes[0].set_xlim(0, N_VOLS - 1)
axes[0].tick_params(axis="both", which="major", labelsize=14)
axes[1].plot(
FULLCURVE_TES,
s0based_fullcurve_signal[:, i_vol],
color="red",
alpha=0.5,
label="$S_0$-Driven Decay Curve",
zorder=0,
)
axes[1].plot(
FULLCURVE_TES,
t2sbased_fullcurve_signal[:, i_vol],
color="blue",
alpha=0.5,
label="$T_{2}^{*}$-Driven Decay Curve",
zorder=1,
)
axes[1].scatter(
MULTIECHO_TES,
s0based_multiecho_signal[:, i_vol],
color="red",
s=150,
alpha=1.,
label="$S_0$-Driven Multi-Echo Signal",
zorder=2,
)
axes[1].scatter(
MULTIECHO_TES,
t2sbased_multiecho_signal[:, i_vol],
color="blue",
s=150,
alpha=1.,
label="$T_{2}^{*}$-Driven Multi-Echo Signal",
zorder=3,
)
axes[1].legend(loc="upper right", fontsize=20)
axes[1].set_ylabel("Signal", fontsize=24)
axes[1].set_xlabel("Echo Time (ms)", fontsize=24)
axes[1].set_xticks(MULTIECHO_TES)
axes[1].set_ylim(0, np.ceil(np.max(s0based_fullcurve_signal) / 1000) * 1000)
axes[1].set_xlim(0, np.max(FULLCURVE_TES))
axes[1].tick_params(axis="both", which="major", labelsize=14)
fig.tight_layout()
# save frame
fig.savefig(filename)
plt.close(fig)
filenames.append(filename)
# build gif
with imageio.get_writer(out_file, mode="I") as writer:
for filename in filenames:
image = imageio.imread(filename)
writer.append_data(image)
# Remove files
for filename in filenames:
os.remove(filename)
with open(out_file, "rb") as file:
display(Image(file.read(), width=600))
Plot \(T_{2}^{*}\) against BOLD signal from single-echo data (TE=30ms)¶
When the BOLD signal is driven entirely by T2* fluctuations, the BOLD signal is very similar to a scaled version of those T2* fluctuations, but there are small differences.
fig, ax = plt.subplots(figsize=(16, 4))
scalar = np.linalg.lstsq(
t2s_ts[:, None],
t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
rcond=None
)[0]
ax.plot(
t2sbased_fullcurve_signal[SINGLEECHO_TE[0], :],
color="orange",
label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
t2s_ts * scalar,
label="$T_{2}^{*}$ (scaled to signal)",
linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
fig.show()
Plot \(S_{0}\) against BOLD signal from single-echo data (TE=30ms)¶
When the BOLD signal is driven entirely by S0 fluctuations, the BOLD signal ends up being a scaled version of the S0 fluctuations.
fig, ax = plt.subplots(figsize=(16, 4))
scalar = np.linalg.lstsq(
s0_ts[:, None],
s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
rcond=None
)[0]
ax.plot(
s0based_fullcurve_signal[SINGLEECHO_TE[0], :],
color="orange",
label=f"BOLD signal (TE={SINGLEECHO_TE[0]}ms)",
)
ax.plot(
s0_ts * scalar,
label="S0 (scaled to signal)",
linestyle="--",
)
ax.set_xlim(0, N_VOLS - 1)
leg = ax.legend()
fig.show()